//    Copyright (C) Kherty.  All rights reserved.
using System;
using System.Diagnostics.Contracts;

namespace OpenLS.Spreadsheet.StandardAddIn.Statistics
{
    class BrentSolver: Solver
    {
        public BrentSolver(DoubleFunction f): base(f, 10000, 1e-6)
        {

        }

      
        public override double Solve(double a, double b)
        {
            double fa = f(a);
            double fb = f(b);
            if (fa * fb > 0)
                throw new ArgumentOutOfRangeException("f(a) * f(b) > 0");
            if (Math.Abs(fa) < Math.Abs(fb))
                return Solve(b, a);
            Contract.Assert(Math.Abs(fa) >= Math.Abs(fb));
            double c = a;
            bool mflag = true;
            double d = b;
            double delta = AbsoluteAccuracy;
            int maximalIterationCount = this.MaximalIterationCount;
            double absoluteAccuracy = this.AbsoluteAccuracy;
            double functionValueAccuracy = this.FunctionValueAccuracy;
            for (int i = 0; i < maximalIterationCount; i++)
            {
                if (Math.Abs(b - a) < absoluteAccuracy)
                    return b;
                if (Math.Abs(fb) < functionValueAccuracy)
                    return b;
                double fc = f(c);
                double s;
                if (fa != fb && fa != fc)
                    s = a * fa * fc / ((fa - fb) * (fa - fc)) + b * fa * fc / ((fb - fa) * (fb - fc)) + c * fa * fb / ((fc - fa) * (fc - fb));
                else
                {
                    s = b - fb * (b - a) / (fb - fa);
                }
                bool condition1 = !isBetween(s, (3 * a + b) / 4, b);
                bool condition2 = mflag && Math.Abs(s - b) >= Math.Abs(b - c) / 2;
                bool condition3 = !mflag && Math.Abs(s - b) >= Math.Abs(c - d) / 2; //\\!!d?
                bool condition4 = mflag && Math.Abs(b - c) < delta;
                bool condition5 = !mflag && Math.Abs(c - d) < delta;
                if (condition1 || condition2 || condition3 || condition4 || condition5)
                {
                    s = (a + b) / 2;
                    mflag = true;
                }
                else
                {
                    mflag = false;
                }
                double fs = f(s);
                d = c;
                if (fa * fs < 0)
                {
                    b = s;
                    fb = f(b);
                }
                else
                {
                    a = s;
                    fa = f(a);
                }
                if (Math.Abs(fa) < Math.Abs(fb))
                {
                    double t = a;
                    a = b;
                    b = t;
                    t = fa;
                    fa = fb;
                    fb = t;
                }


            }
            throw new InvalidOperationException("Convergence");
        }

        private bool isBetween(double number, double bound1, double bound2)
        {
            return (number - bound1) * (number - bound2) <= 0;
        }

    }
}